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Abstract 

A  model  for  planar  solid  oxide  fuel  cell  repeat  elements  and  stacks  has  been  developed.  Distribution  of  concentrations,  reaction  rates 
and  temperatures  (both  gases  and  solids)  are  computed  as  well  as  overall  performance  results.  Specific  experiments  provide  inputs  to  the 
model  by  a  parameter  estimation  method. 

The  modeling  approach  developed  allows  to  compare  several  configurations.  As  the  number  of  design  parameters  is  large  (from  cell 
size,  component  thicknesses  to  gas  flow  configuration),  the  model  is  designed  to  change  easily  these  parameters  so  as  to  explore  as  many 
cases  as  possible.  This  is  particularly  true  for  the  flow  configuration  (inlet  position,  outlets)  for  which  several  options  are  considered. 

This  model  assists  in  choosing  a  configuration  and  allows  to  perform  sensitivity  studies  in  an  efficient  way  (without  having  to  produce  a 
new  mesh  such  as  for  CFD  tools)  or  to  be  combined  with  an  optimization  tool.  A  first  validation  with  experimental  results,  performed  on 
a  particular  stack  design,  is  presented.  Issues  of  model  accuracy  and  sensitivity  to  uncertain  inputs  are  discussed. 

©  2004  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Design  and  development  of  a  solid  oxide  fuel  cell  stack 
is  a  task  involving  many  different  aspects.  The  goals  for 
a  stack  can  be  very  challenging:  maximum  power  output, 
compactness  and  long  term  stability.  The  trend  to  increase 
power  densities  may  be  in  contradiction  with  degradation 
behavior,  the  drawback  being  an  increased  temperature  and 
gradient,  which  is  likely  to  induce  accelerated  degradation 
and  failure.  On  the  other  side,  constraints  from  the  system 
(such  as  the  maximum  pressure  drop  on  the  stack),  ceramic 
fabrication  feasibility  and  cost  aspects  add  more  elements 
to  consider. 

Modeling  can  be  a  tool  to  help  decision  making  on  some 
important  stack  characteristics.  Several  models  are  published 
in  literature  [1-3].  These  models  represent  a  particular  de¬ 
fined  configuration. 

This  work  is  being  realized  in  the  frame  of  a  stack  devel¬ 
opment  project.  The  focus  and  aim  of  the  model  developed 
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is  to  be  able  to  compare  different  configuration,  explore 
geometrical  parameters  and  predict  performance  and  behav¬ 
ior  of  the  existing  repeat  element  and  stacks. 

The  complete  mathematical  description  is  based  on  the 
following  equations:  conservation  of  momentum  with  a 
porous  media  description,  conservation  of  species  for  both 
gases  (expressed  in  moles),  conservation  of  energy  for 
gases,  conservation  of  energy  for  the  solid  (assumed  as 
monolithic). 

This  work  assumes  that  a  2D  model,  combined  with 
optimization,  will  allow  to  identify  better  configurations, 
optimize  the  decision  variables  and  reduce  the  number  of 
experiments  in  the  development  process.  In  the  different 
levels  of  details  for  stack  modeling,  CFD  models  are  accu¬ 
rate  with  a  good  geometry  definition  and  complete  solving 
of  transport  phenomena,  on  the  other  hand,  a  ID  model 
would  not  allow  to  capture  the  behavior  of  the  considered 
configuration.  The  gap  between  these  two  model  types  has 
to  be  filled  in,  as  CFD  is  not  an  efficient  tool  to  explore 
new  geometry  and  optimize  design  parameters:  the  mesh 
generation  is  a  long  process  and  CPU  time  for  simulation 
is  long.  The  presented  model  aims  to  provide  the  essential 
information  from  a  CFD  model  (such  as  velocity,  pressure, 
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concentration  and  temperature  field)  with  a  more  efficient 
computing  time  and  flexibility. 

Equations  and  assumptions  made  in  the  model  are  pre¬ 
sented,  then  important  inputs  and  results  obtained  are  shown. 
Accuracy  issues  are  discussed.  As  the  model  needs  to  be 
validated  experimentally  to  be  used  as  a  decision  tool,  some 
elements  of  validation  are  presented.  Finally,  an  example  of 
comparison  of  configurations  is  presented. 

2.  The  repeat  element  configuration 

The  repeat  element  considered  is  based  on  an  anode  sup¬ 
ported  cell  and  metallic  interconnects  and  is  designed  to 
operate  around  1050  K.  This  concept  has  demonstrated  an 
electrical  output,  100  We  with  a  stack  of  six  cells  (of  52  cm2 
active  area  per  cell).  One  of  the  main  features  of  this  stack 
is  to  used  internal  manifolding  feed  the  reactants  to  the  cell, 
therefore  no  additional  pieces  are  needed  and  the  stack  is 
very  compact.  The  inlet  of  gases  is  punctual  and  the  outlet  is 
distributed  on  a  side.  The  current  lines  of  the  flow  field  are 
therefore  in  2D  (see  in  Fig.  la).  On  most  of  the  cell  surface, 
the  flow  is  approximately  in  counter  flow  mode. 

The  current  repeat  element  uses  square  cells  of  80  mm 
side,  the  active  surface  (accounting  for  the  surface  occupied 
by  sealing  and  border)  is  of  ca.  52  cm2.  The  total  repeat 
element  thickness  is  below  3  mm  with  a  250  jam  thick  cell 
(200  [xm  for  the  anode  support,  6  jam  thin  electrolyte  and  ca. 


50  |xm  porous  cathode),  interconnects  of  0.75  mm  and  gas 
diffusion  layers  of  0.5  mm  on  both  sides. 

From  the  experimental  side,  the  different  dimensions  are 
fixed.  In  a  design  mode,  however,  the  degrees  of  freedom 
are:  the  cell  total  area  and  geometry  (length,  aspect  ratio), 
the  inlet  position,  the  gas  diffusion  layer  and  interconnect 
thickness.  These  geometrical  decision  variables  are  linked 
to  the  operating  variables  such  as  the  fuel  and  air  flow  rates 
(under  the  constraint  of  pressure  drop  targets  for  the  repeat 
element). 

Based  on  the  same  technology  using  internal  manifolding, 
other  configurations  are  possible,  one  example  is  a  case  with 
a  central  feed  of  hydrogen  and  two  inlets  for  the  air  flow 
which  represented  in  Fig.  1(b).  A  comparison  of  these  two 
configurations  is  presented  in  Section  6. 

3.  Model 

The  repeat  element  configurations  considered  are  in  ob¬ 
vious  need  for  2D  modeling.  The  fluid  patterns  cannot  be 
represented  by  simple  mono-dimensional  flow.  The  geomet¬ 
ric  definition  of  the  model  is  presented  in  Fig.  1(c)  for  the 
present  configuration.  As  the  cell  design  is  symmetric  (with 
an  axes  passing  through  the  inlets),  only  half  of  the  cell 
have  been  modeled.  The  lengths  of  the  cell  are  variables 
of  the  model,  the  inlet  coordinates  are  parameters  of  the 
model. 
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Fig.  1.  Configuration  simulated  and  geometrical  models,  (a)  Configuration  “counter  flow,”  (b)  configuration  central  feed,  (c)  model  for  “counter  flow” 
geometry,  (d)  model  for  central  feed  geometry. 


306 


D.  Larrain  et  al.  /  Journal  of  Power  Sources  131  (2004)  304-312 


In  Fig.  1(d),  two  axes  of  symmetry  (one  passing  by  the 
three  inlets,  and  the  other  one  by  the  central  inlet)  allow 
to  model  only  a  quarter  of  the  cell.  The  following  section 
presents  the  assumptions  and  equations  used  in  the  model. 

3.1.  Fluid  flow  and  species  balance 

The  fluid  pattern  in  the  gas  diffusion  layer  is  a  3D  phe¬ 
nomenon  with  2D  streamlines  and  a  velocity  profile  in  the 
height  of  the  channel.  The  model  developed  assumes  a  plug 
flow  in  the  height  of  the  channel  to  reduce  the  flow  descrip¬ 
tion  to  a  2D  flow  pattern,  which  is  the  main  feature  of  this 
repeat  element  model. 

To  describe  the  fluid  flow  and  species  balance  equa¬ 
tions,  several  options  are  possible.  The  main  assumptions 
made  here  are:  the  total  molar  concentration  of  species  is 
decoupled  from  the  temperature  field,  the  change  in  to¬ 
tal  number  of  moles  is  decoupled  from  the  velocity  field, 
and  the  viscosity  of  the  fluid  is  assumed  constant  as  well. 
These  three  assumptions  decouple  the  fluid  motion  from  the 
energy  equations  and  the  species  balance.  The  essential 
coupling  of  the  reaction  rate  with  the  species  concentration 
and  temperature  is  made  by  the  kinetic  model  used.  Solving 
the  complete  set  of  equations  for  the  momentum,  species 
balance  and  energy  requires  a  CFD  tool  using  a  finite  vol¬ 
ume  method  for  the  resolution  which  is  not  the  goal  of  the 
present  work. 

If  no  reforming  is  considered  (operation  with  hydrogen) 
the  total  number  of  moles  is  constant  on  the  fuel  side. 
However,  on  the  air  side  the  number  of  species  decreases: 
considering  for  example  an  air  ratio  of  2  and  70%  fuel  uti¬ 
lization,  the  oxygen  utilization  is  35%  and  the  total  number 
of  moles  decreases  by  ca.  8%.  Then  the  error  on  the  veloc¬ 
ity  field  for  the  air  side  is  of  ca.  8%.  The  decoupling  from 
the  temperature  causes  the  same  amount  of  error  to  the 
velocity  field  with  a  maximum  temperature  difference  of 
100  K  (which  is  quite  common  in  SOFC).  Nevertheless,  as 
the  species  balance  is  expressed  in  terms  of  molar  flux,  this 
assumption  does  not  affect  the  consistency  of  the  species 
balance.  The  main  error  expected  is  on  the  pressure  drop 
estimation  which  does  not  account  for  the  acceleration  of 
the  fluid  neither  for  the  change  in  viscosity. 

The  presented  assumptions  lead  to  the  following  equa¬ 
tion  system.  For  the  fluid  flow,  the  velocity  is  linked  to  the 
pressure  field  by  a  Darcy  equation  in  2D. 

-VC P)  =  V  (1) 

where  v  is  the  velocity  vector  on  v  and  y\  p  is  the  dynamic 
viscosity;  P  is  the  pressure;  K  is  the  permeability  coefficient 
determined  experimentally. 

Since  the  total  mole  change  is  assumed  zero  and  the  vis¬ 
cosity  constant,  the  pressure  field  is  computed  from  the 
boundary  conditions  with  a  Laplace  equation 


Finally,  the  species  balance  is  therefore  computed  with  an 
equation  including  convective  and  diffusive  transport: 

V  •  Fi  —  DACj  =  rt  (3) 

where  Fi  is  the  local  molar  flux  vector  (on  v  and  y  direc¬ 
tion)  of  species  /;  f;  is  the  rate  of  reaction  per  unit  volume; 
D  is  the  binary  diffusion  coefficient  computed  from  the 
Fuller-Schettler-Giddings  equation  [4]. 

The  boundary  conditions  applied  on  the  fluid  equations 
are: 

•  inlet  concentration  fixed  at  inlet 

•  uniform  pressure  at  the  outlet 

•  as  the  punctual  inlet  is  a  mathematical  singularity  (infinite 
velocities),  the  rate  of  reaction  and  the  velocity  are  fixed 
to  zero  at  the  inlet 

•  the  wall  boundary  condition  is  defined  setting  the  velocity 
vector  component  normal  to  the  wall  to  zero,  and  the 
diffusive  term  ( dCi/dn )  =0. 

•  as  post-combustion  is  included  in  the  model,  the  domain 
modeling  the  fuel  flow  is  extended  to  a  post-combustion 
area  where  an  almost  complete  combustion  of  hydrogen 
is  assumed  at  the  outlet  of  the  domain. 


3.2.  Energy  equations 


Energy  equations  are  defined  for  the  solid  part,  assuming 
averaged  thermal  properties  [1,2, 5, 6].  The  method  used  for 
volume  averaging  here  is  a  simple  connexion  of  series  and 
parallel  conduction  [7,8],  the  thermal  transport  properties 
are  considered  as  isotropic.  The  solid  energy  equation  is 
therefore: 


+  2  =  0 


where  kSx  is  the  average  thermal  conductivity,  0  is  the 
sum  of  the  volumetric  sources  including  the  heat  from  the 
electrochemical  reaction,  the  electric  power  removed  from 
the  system  and  heat  transfer  to  the  fluids.  For  the  fluids,  the 
energy  equation  is  a  local  enthalpy  balance: 
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where  CPi  is  the  heat  capacity  of  species  /;  Fi-X  is  the  local 
molar  flux  of  i  in  the  v  direction;  and  0  is  the  volumetric 
heat  sources  which  are  essentially  the  heat  transfer  with  the 
solid  and  the  enthalpy  of  reactants  and  products. 


3.3.  Kinetics 


As  the  model  represents  a  repeat  element  based  on  anode 
supported  cells,  the  current  path  through  the  electrolyte  is 
assumed  to  be  normal  to  the  surface  and  the  equation  for 
the  conservation  of  charge  is  not  solved.  The  local  reaction 
rate  is  computed  with  the  following  scheme: 

f^cell  =  f^Nernst  ^ohmic  x  j  77 act 


AP  =  0 


(6) 
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where  t/Nernst  is  the  Nernst  potential  computed  locally  from 
the  temperature  and  the  partial  pressures;  UCQ n  is  the  cell 
operating  voltage;  j  is  the  local  current  density;  and  r/aa  is 
the  sum  of  overpotentials.  The  local  resistance  and  activa¬ 
tion  overpotential  are  functions  of  temperature  and  current 
density. 

3.4.  Implementation  of  the  model 

The  set  of  equation  described  has  been  implemented 
in  the  gPROMS  tool  which  is  an  equation  solver  based 
tool  allowing  to  solve  distributed  domains  [9].  The 
scheme  used  to  solve  the  distributed  domain  is  cen¬ 
tered  finite  difference.  This  tool  have  embedded  param¬ 
eter  estimation  and  optimization  solvers.  As  this  tool  is 
based  on  a  equation  solver,  it  gives  the  possibility  to 
change  the  input  variables  and  degrees  of  freedom  of  the 
problem. 

The  presented  equations  have  been  implemented  in  a  nor¬ 
malized  form  allowing  efficient  sensitivity  on  the  shape  and 
size  of  the  cell. 


4.  Model  outputs  and  validation 

The  results  from  the  implemented  model  are  dependent  on 
the  boundary  conditions  applied  and  the  input  variables  used 
for  some  properties.  The  outputs  of  the  model  are  presented 
and  validation  of  the  model  consistency  with  the  mesh  size 
is  discussed. 

4.1.  Thermal  boundary  conditions  applied 

Boundary  conditions  applied  to  the  model  are  essential 
as  they  affect  the  temperature  field.  Usually  repeat  element 
are  assumed  to  be  in  the  middle  of  a  stack  with  adiabatic 
boundary  condition.  This  boundary  condition  is  easily  de¬ 
fined.  However,  previous  work  [1,10]  shows  that  for  a  stack 
mounted  with  metallic  interconnects,  no  adiabatic  cells  are 
found.  Adiabatic  conditions  are  then  the  worst  case  scenario, 
but  do  not  represent  the  reality  for  short  stacks.  The  results 
for  the  IV  characteristics  presented  are  simulated  with  heat 
flux  to  the  environment  to  represent  a  single  repeat  element 
test. 

Post-combustion  is  considered  in  the  model.  To  com¬ 
pute  this  contribution,  the  enthalpy  of  combustion  of  the 
remaining  species  of  hydrogen  is  computed  on  each  point 
of  the  fuel  outlet  (see  in  Fig.  1(c))  [2].  However,  one  of 
the  uncertainties  on  the  post-combustion  contribution  is 
the  proportion  of  heat  which  is  really  absorbed  by  the  re¬ 
peat  element,  this  depends  on  where  the  post-combustion 
flame  is  located.  This  flame  starts  on  the  border  of  the 
stack  but  extends  to  the  environment.  Therefore,  it  is 
assumed  that  part  of  the  heat  only  is  absorbed  by  the 
stack. 


4.2.  Model  inputs 

The  equations  defined  in  Section  3  need  some  input 
variables  to  be  defined.  When  available,  experimental  data 
is  used,  however  some  values  are  known  and  found  in  the 
literature.  Some  of  these  inputs  are  uncertain  and  the  re¬ 
sponse  has  to  be  evaluated. 

The  heat  conduction  of  the  cell  is  computed  from  the 
value  for  the  anode  cermet  [11],  electrolyte  and  cathode. 
The  cermet  is  assumed  to  have  40%  volume  Ni  and  50% 
porosity.  Recent  values  for  perovskites  materials  (cathode) 
are  difficult  to  find  [12]. 

The  heat  transfer  coefficient  is  computed  from  a  Nusselt 
number  for  laminar  flow  between  parallel  plates  (^  8)  and 
the  exchange  surface  takes  into  account  the  surface  added 
by  the  current  collector.  Some  uncertainties  remain  on  the 
heat  transfer  coefficient,  especially  when  considering  that 
the  surface  layers  of  the  cell  are  porous,  but  the  model  is  not 
very  sensitive  to  this  parameter  (for  Nusselt  numbers  in  the 
range  of  8-20,  the  magnitude  of  the  temperature  response 
differs  by  10  K  at  the  maximum  temperature  point). 

Another  uncertainty  is  the  proportion  of  heat  from  the 
post-combustion  which  is  absorbed  by  the  stack.  However, 
it  has  a  small  incidence  on  the  performance  and  temperature 
results  if  varied  in  the  range  0.25-0.75. 

For  the  radiative  boundary  condition,  the  emissivity  is 
assumed  to  be  at  0.9.  But  again,  varying  this  parameter 
between  0.7  and  0.9  has  not  a  major  influence  on  the  results. 

The  input  values  for  the  kinetic  parameters  are  obtained  by 
a  parameter  identification  method  from  experimental  results 
on  a  button  cell  of  1  cm2  active  area  [13].  The  contribution 
of  the  interface  resistance  between  current  collectors  and 
interconnects  are  taken  from  experiments  [14]  (Table  1)  . 

The  lengths  of  the  cell  edges  are  variables  and  therefore 
sensitivity  to  the  cell  size  and  to  the  aspect  ratio  can  easily 
be  performed. 

The  decision  variables  have  been  listed  in  Section  2. 

4.3.  Model  outputs 

The  model  computes  fields  of  pressure,  velocity,  con¬ 
centration,  reaction  rates  and  temperature.  Fig.  2  shows 
the  concentration  field  for  hydrogen  and  solid  temperature. 
From  Fig.  2(b),  it  can  be  noticed  that  punctual  inlet  induces 
large  regions  with  lean  fuel  concentration.  At  the  outlet, 
the  concentration  profile  is  not  homogeneous,  with  higher 

Table  1 


Input  variables  for  the  model 


Input 

Units 

Range 

Cermet  conductivity 

W/mK 

8.9 

Interconnect  conductivity 

W/mK 

25  (@1100  K) 

Fraction  heat  post-combustion 

— 

0.25-0.75  (0.5) 

Emissivity 

— 

0.9 

Nusselt  number  (two  plates) 

— 

O 

~  O 
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Fig.  2.  Model  outputs,  all  cases  @28  A  total  current  and  70%  fuel  utilization, 
mole  fraction. 

concentration  close  to  the  middle  of  the  cell.  The  steep  de¬ 
crease  in  concentration  near  the  inlet  is  due  to  the  extremely 
high  current  density  at  the  inlet. 

For  the  temperature  field  (seen  in  Fig.  2(a)),  the  hot  spot 
is  located  close  to  the  fuel  entrance  as  expected  in  a  counter 
flow  situation.  As  the  air  flux  is  low,  the  hot  spot  is  between 
the  two  inlets  and  the  general  temperature  profile  is  sim¬ 
ilar  for  adiabatic  and  non-adiabatic  boundary  conditions. 
The  maximum  temperature  changes  from  1240  to  1085  K 
when  changing  the  boundary  conditions  from  adiabatic  to 
non-adiabatic. 

4.4.  Validity  of  the  model  and  accuracy  problems 

The  results  presented  are  computed  with  a  finite  difference 
scheme  and  the  accuracy  of  the  results  can  be  affected  by  the 
size  of  the  mesh  used.  This  is  specially  the  case  here  where 
a  rather  coarse  mesh  is  currently  used  (with  200-1500  nodes 
on  the  surface). 

Fig.  3(a)  shows  how  the  results  on  the  temperature  profile 
can  be  affected.  The  mesh  size  was  varied  from  11  x  21  to 


Y 


(a)  Temperature  solid  adiabatic,  maximum  temperature  1240  K,  (b)  hydrogen 

31  x  61  to  study  the  sensitivity  of  the  results.  As  expected,  the 
coarser  mesh  has  a  maximum  temperature  higher  than  the 
finer  one  (of  ca.  25  K).  Nevertheless,  the  shape  of  the  tem¬ 
perature  profile  remains  the  same.  It  can  be  pointed  out  that 
for  the  two  finer  meshes,  the  temperature  profile  is  similar. 

The  error  on  the  species  balance  decreases  with  increasing 
mesh  size  from  1.2%  on  the  coarser  to  less  than  0.3%  as 
reported  in  Table  2.  The  error  seems  to  be  related  to  the 
mesh  definition.  As  the  change  in  intensity  and  direction  of 
velocity  are  very  steep  close  to  the  inlet,  most  of  the  error 
may  come  from  the  inlet  region.  Ideally,  the  mesh  should  be 
refined  near  the  inlet  but  this  is  unfortunately  not  possible 
with  the  tool  used. 

In  order  to  use  the  model  for  optimization  of  decision 
variables,  the  important  point  is  to  verify  that  the  mesh  size 
has  no  impact  on  the  sensitivity  of  the  chosen  criteria.  For 
optimization,  two  objective  functions  are  considered:  maxi¬ 
mizing  the  specific  power  output  per  unit  volume  in  W/cm3 
and  minimizing  the  maximum  temperature  of  the  solid.  In 
Fig.  3(b),  the  sensitivity  to  the  interconnect  thickness  of 
those  two  criteria  is  shown  for  two  different  mesh  sizes.  The 
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Fig.  3.  Sensitivity  to  mesh  size:  temperature  profile  and  response  to  a  parameter  change,  (a)  Temperature  profile  variation  with  increasing  mesh  size, 
@28  A  total  current,  70%  f.u.,  adiabatic,  the  CPU  time  was  2,  15,  30  and  600  min.  (b)  Sensitivity  to  interconnect  thickness  for  two  different  mesh  sizes 
(mesh  1:  11  x  21,  mesh  2:  16  x  31). 


Table  2 

Results  sensitivity  to  the  mesh  size 
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11  x  21 

16  x  31 

21  x  41 

31  x  61 

41  x  81 

51  x  101 

Error  on  species  balance  (%) 

1.2 

0.6 

0.45 

0.38 

0.32 

0.28 

Maximum  temperature  solid  (K) 

1240.2 

1224.5 

1217.2 

1216.6 

— 

— 

Power  output  (W) 

19.80 

19.92 

19.92 

19.93 

— 

— 

absolute  value  of  the  temperature  differs  by  ca.  15  K  but 
the  response  is  perfectly  similar.  The  model  developed  can 
therefore  be  used  for  optimization  even  with  a  coarse  grid 
(which  is  required  in  term  of  CPU  time). 

5.  Experimental  validation 

The  developed  model  has  been  compared  to  experimental 
results  for  a  repeat  element.  The  performance  is  compared 
for  three  different  operating  conditions  (in  terms  of  hydro¬ 
gen  flux).  The  operating  temperature  was  1045  K  (i.e.  the 
environment  temperature  for  the  model)  measured  in  the  test 
oven. 

5.7.  IV  characteristic  validation 

The  results  presented  in  Fig.  4  show  simulated  and  mea¬ 
sured  IV  curves  for  the  three  cases  (@200,  300,  360ml/min 
hydrogen  flux  and  air  ratio  of  2.5).  The  cell  was  run  up  to 
68-70%  fuel  utilization.  The  simulated  data  uses  kinetic 
parameters  identified  on  button  cells  (see  Section  3.3).  The 
model  should  then  be  able  to  simulate  the  performance  of 
a  repeat  element  using  the  same  cells.  The  comparison  of 
simulated  and  experimental  characteristics  shows  that  the 
response  of  the  model  to  the  flux  is  satisfactory,  but  that  there 
is  a  large  error  on  the  open  circuit  voltage.  The  simulated 
open  circuit  voltage  (OCV)  is  overestimated  by  ca.  100  mV. 

The  conclusion  here  is  that  the  model  used  for  simula¬ 
tion  does  not  take  into  account  a  phenomena  responsible  for 


Fig.  4.  Current  potential  (IV)  characteristics:  simulated  and  experimental; 
52  cm2,  environment  @1045K,  air  ratio  2.5. 


lower  OCV.  This  lower  OCV  can  be  attributed  to  several 
phenomena:  diffusion  of  species  from  the  post-combustion 
area,  leakage  current  [15]  and  problems  on  the  sealing  cre¬ 
ating  a  cross-over  of  reactants. 

In  previous  work,  implementing  diffusion  from  the 
post-combustion  zone  in  a  model  for  a  cell  based  on  elec¬ 
trolyte  supported  cells,  good  simulation  of  OCV  depending 
on  flux  was  obtained  [2] .  This  aspect  is  implemented  in  the 
model  (see  Section  3.1)  and  contributes  to  a  lowering  of 
the  OCV  in  the  order  of  20-40  mV  compared  to  theoretical 
Nernst  voltage.  Leakage  current  has  therefore  to  be  added 
to  the  model. 

5.2.  Leakage  current  model 

At  OCV,  assuming  the  electrolyte  is  a  pure  ionic  conduc¬ 
tor,  no  consumption  of  species  neither  ionic  current  should 
be  found.  But  if  the  electrolyte  has  some  electronic  conduc¬ 
tivity,  then  it  behaves  like  a  mixed  conductor  and  some  ionic 
current  exists  even  at  OCV.  When  the  cell  is  supplying  cur¬ 
rent  to  a  load,  the  external  current  is  added  to  the  internal 
short  circuit  current  through  the  electrolyte.  In  anode  sup¬ 
ported  cells,  the  electrolyte  is  very  thin  (6  pirn  in  our  case). 
From  literature,  the  electronic  resistance  for  this  thickness 
is  in  the  range  of  5-100  Q  cm2  [16-18]. 

The  model  chosen  to  represent  the  leakage  current  is  based 
on  the  equivalent  circuit  shown  in  Fig.  5.  At  OCV  the  leak¬ 
age  current  is  essentially  determined  from  the  ionic  and  elec¬ 
tronic  resistance.  The  ionic  resistance  for  the  cell  used  in 
our  experiment  is  assumed  to  be  known  [15],  with  a  value 
of  0.12  £2  cm2. 

The  equations  describing  the  system  are: 

E0  ~  ^ionicjion  —  ^totjload  —  Ucen  (7) 

Jload  T  Jloss  —  Jion  (8) 


Eo 


Fig.  5.  Equivalent  circuit  for  the  electrode/electrolyte  system. 
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Table  3 


Cases  for  the  OCV  simulation  comparison 


Fluxes 

Centered 

Counter 

Cel  Ce2 

Ce3 

Col 

Co2 

Co3 

Fuel  flow  (ml/min) 

180  230 

260 

250 

300 

400 

Air  flow  (1/min) 

2  2 

2 

2.3 

2.5 

2.5 

where  E 0  is  the  local  Nernst  voltage;  Tronic  is  the  ionic 
resistance  of  the  electrolyte;  jion  is  the  ionic  current  through 
the  electrolyte;  Rtoi  is  the  resistance  of  the  electrodes  and 
interconnects;  jioad  is  the  current  density  in  the  electrodes; 
UCq it  is  the  cell  operating  voltage. 

At  OCV,  the  current  in  the  external  circuit  is  zero  and 
therefore  the  loss  current  can  be  computed  at  OCV  as: 


^ionic  +  ^elec 

where  Re \QC  is  the  electronic  resistance  of  the  electrolyte. 
For  any  operating  point,  the  local  loss  current  is  described 
by  the  equation: 

^elecjloss  —  f^cell  T  ^tot  Jioad  (10) 

The  value  of  the  electronic  resistance  has  to  be  estimated. 
Parameter  identification  method  has  been  used  to  find  the 
estimate.  The  experimental  data  used  for  the  parameter  es¬ 
timation  comes  from  three  different  experiments: 

•  a  button  cell  of  16  cm2  area  with  central  feed  [13] 

•  a  repeat  element  with  central  feed  of  air  and  fuel  on  a 
square  cell  (80  mm  side) 

•  a  regular  repeat  element  with  improved  sealing  (compared 
to  the  IV  curve  experiment). 

For  these  three  experiments,  three  sets  of  conditions  were 
available  (see  in  Table  3).  The  identification  has  been  per¬ 
formed  on  the  three  experiments  with  a  common  assumption 
that  the  electrolyte  thickness  was  of  6  [xm  (the  thickness  has 
not  been  measured).  The  result  from  the  parameter  estima¬ 
tion  is  a  value  of  9.58  £2  cm2  for  the  electrolyte  electronic 


resistance.  The  95%  confidence  interval  is  quite  large  with  a 
range  of  ±2  Q  cm2,  the  poor  statistical  quality  of  the  result 
is  due  to  a  lack  of  data. 

The  comparison  of  the  OCV  on  the  two  different  repeat 
element  configurations  is  presented  in  Fig.  6.  The  application 
of  the  leakage  current  model  and  the  parameter  optimization 
allow  a  sensible  improvement  of  the  OCV  simulation.  The 
remaining  error  is  below  20  mV.  The  quality  of  the  model 
could  be  improved  with  further  experiments  and  parameter 
identification. 

5.3.  New  simulation  of  the  IV  curve 

The  model  for  OCV  calculation  has  then  been  applied  to 
simulate  the  IV  characteristics.  With  the  identified  value,  an 
offset  of  ca.  40-50  mV  remains.  This  could  be  explained  by 
the  fact  that  some  reactants  cross-over  occurs  in  this  repeat 
element  (as  the  experimental  OCV  does  not  increase  with  the 
flux  in  this  case)  or  to  a  slightly  lower  electrolyte  thickness. 
Correcting  the  electrolyte  thickness  to  4.5  ptm  in  the  model 
inputs,  the  remaining  error  is  of  ca.  20  mV  as  shown  in 
Fig.  7,  which  is  satisfactory. 

This  last  result  shows  that  comparison  of  experimental 
data  with  simulation  is  essential  to  improve  models.  In  this 
case,  the  error  in  OCV  has  necessitated  to  add  a  new  phe¬ 
nomenon  into  the  model  and  the  order  of  magnitude  of  this 
phenomenon  has  been  identified. 

From  the  identified  values  for  the  electronic  conductivity, 
the  leakage  current  is  in  the  order  of  0. 1  A/cm2  which  is 
not  negligible  at  all.  This  affects  the  species  balance  as 
the  leakage  current  involves  conversion  of  hydrogen,  there¬ 
fore  the  effective  fuel  utilization  is  higher  than  the  fuel 
utilization  computed  from  the  current  drawn  in  the  exter¬ 
nal  circuit.  This  is  particularly  important  and  critical  at 
high  fuel  utilization.  When  operating  at  70%  fuel  utiliza¬ 
tion,  the  conditions  in  the  cell  are  equivalent  to  operation 
at  75  or  80%.  The  leakage  current  may  therefore  limit  the 
efficiency  of  the  stack  by  limiting  operation  at  high  fuel 
utilization. 


Fig.  6.  Experimental  vs.  simulated  OCV  with  estimated  value  for  Re iec,  filled  sign  are  the  estimated  value,  stars  and  cross  gives  simulation  results  with 
the  lower  and  higher  bound  of  the  95%  confidence  interval,  (a)  Central  feed  repeat  element,  (b)  “counter”  flow  repeat  element. 
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Fig.  7.  Current  potential  (IV)  characteristics:  simulated  (with  leakage 
current  model)  and  experimental;  52  cm2,  environment  @1045K,  air  ra¬ 
tio  2.5. 

6.  Comparison  of  configurations 

The  counter  flow  configuration  (in  Fig.  1(c))  is  the  first 
one  to  be  experimentally  demonstrated  so  far,  but  other  con¬ 
figurations  are  possible  based  on  the  internal  manifolding 
concept.  The  example  shown  next  is  described  in  Fig.  1(d). 
It  uses  a  central  feed  for  the  fuel  and  two  inlets  for  the  air 
flow,  the  flow  pattern  being  closer  to  a  co-flow.  To  com¬ 
pare  the  results  for  this  configuration  with  the  counter  flow 
case,  simulations  have  been  performed  with  the  same  cell 
size,  same  thickness  for  the  interconnect,  the  same  channel 
height  and  operating  conditions  (fluxes  and  temperature). 
The  electrochemical  characteristic  is  compared  in  parallel  to 
the  maximum  temperature. 

Fig.  8  shows  the  electrochemical  performance  to  be  quite 
similar  for  both  cases.  The  maximum  temperature  encoun¬ 
tered  in  the  solid  is  quite  different.  The  temperature  is  much 
higher  in  the  counter  flow  configuration.  This  can  be  ex¬ 
plained  by  the  post-combustion  that  occurs  only  on  one 
side  for  the  counter  flow  case.  The  intensity  of  the  heat 
release  is  then  much  higher  than  in  the  other  case  as  the 


Fig.  8.  Configuration  comparison  between  “counter”  flow  and  central  feed. 


post-combustion  takes  place  on  the  whole  perimeter  of  the 
cell. 

A  change  in  configuration  can  thus  be  an  interesting  so¬ 
lution  to  achieve  a  similar  power  output  with  lower  tem¬ 
perature  levels.  A  limit  to  the  actual  description  of  the  flow 
field  can  be  found  with  this  case.  To  be  effectively  realized, 
sealing  rings  will  have  to  be  introduced,  changing  the  flow 
pattern  in  the  fuel  side,  this  change  is  here  significant  as  the 
fluid  path  is  short.  This  may  cause  an  overestimation  of  the 
electrochemical  performance  of  such  a  cell. 

7.  Conclusion 

The  developed  modeling  tool  has  demonstrated  ability  to 
predict  experimental  results  for  the  electrochemical  perfor¬ 
mance  of  technical  repeat  element,  to  carry  out  sensitivity 
studies,  and  to  be  adapted  to  simulate  new  configurations. 
Therefore,  in  the  present  status  this  tool  can  be  used  to  op¬ 
timize  repeat  element  and  explore  new  configurations. 

Experimental  validation  of  the  temperature  field  and  on 
further  experiments  is  expected  to  enhance  the  model  and 
confidence  on  the  results.  Comparison  of  the  results  for  the 
same  case  in  a  CFD  model  will  be  carried  out.  CFD  may  be 
used  as  well  to  verify  the  range  of  validity  of  some  assump¬ 
tions. 

Boundary  conditions  definition  needs  to  be  further 
studied:  adiabatic  boundary  conditions  are  the  worst  case 
scenario,  not  expected  to  be  found  in  stack  with  metallic 
interconnects  of  less  than  30  cells.  This  will  have  to  be  veri¬ 
fied  by  using  this  repeat  element  model  as  a  base  for  a  stack 
model.  Then,  the  results  from  this  stack  model  will  have 
to  be  used  to  define  more  realistic  boundary  conditions. 
The  model  will  be  extended  to  hydrocarbon  reformed  fuels 
which  are  planned  be  tested  soon  experimentally.  Finally, 
further  experiments  are  necessary  to  verify  the  preliminary 
results  obtained  for  the  leakage  current. 

As  the  simulation  of  an  operating  point  can  be  quite  fast, 
this  model  will  be  coupled  to  a  multi-objective  optimization 
algorithm  to  explore  the  range  of  possible  solutions  for  each 
configuration  and  identify  the  interesting  solutions.  The 
tool  presented  here  will  be  used  in  parallel  to  other  models 
with  different  levels  of  details  (with  CFD)  and  different 
scales  (stack  model  and  system  model)  to  form  a  conception 
and  design  platform  for  SOFC  planar  repeat  element  and 
stacks. 
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